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We report on conditional statistics of turbulent entrainment in a temporal plane jet at 
Re = 5000 along enstrophy isosurfaces. The data is obtained by direct numerical simula- 
tion and the statistics span 24 orders of magnitude in the enstrophy threshold, ranging 
^ . from essentially irrotational fluid outside the jet to the turbulent core region at the jet 

M-H ' centre. We use two independent estimators for the local entrainment velocity Vn based 

C/3 . on the enstrophy budget and develop a theoretical model that describes the expected be- 

haviour inside the viscous superlayer. The model is in reasonably good agreement with 
£^ . the data and predicts that viscous transport has magnitude 2vn and viscous dissipation 

^H I —Vn- Although Vn IS similar in magnitude to the Kolmogorov velocity m^, its value is not 

Oj. independent of the threshold. Contrarily, the product between Vn and surface area S of 

the viscous superlayer, the entrainment flux, is independent of threshold implying that 
the local transport velocity and surface area adjust concertedly to the imposed global 
<^ ' rate at the fringes of the turbulent layer. We divide the flow into a (i) viscous superlayer 

\j^ . in which the transmission of enstrophy is dominated by viscous effects; (ii) a turbulent 

t^^ ' core region which is dominated by inertial production, viscous transport and dissipation; 

\l . and (iii) an intermediate region that connects the two. The average distance between 

the core region and the viscous superlayer is about 10 Kolmogorov length scales, indi- 
^^ , eating that intense turbulent flow regions and viscosity-dominated regions are in close 

O 

m 



proximity. 



1. Introduction 

• 1-^ 

j^ I Turbulent entrainment is the incorporation of ambient fluid at the boundary of turbu- 

H ' lent flows such as free shear flows or at the free stream edge of turbulent boundary layers. 

- - -' It is an important process in a variety of engineering and geophysical flows controlling the 

turbulent transfer of mass, heat and momentum (Pope 2000; Stull 1998; Thorpe 2005). 
A relevant, yet unresolved issue that has received renewed interest in recent years is the 
connection between processes that are associated with the large scale organization of the 
flow and processes that occur at the scale of the smallest eddies (e.g. Westerweel et al. 
2005; Da Silva & Taveira 2010; Hunt et al. 2011; Philip & Marusic 2012; Wolf et al. 
2012). 

The integral rate at which ambient fluid is incorporated into the turbulent flow is 
independent of the small scale details of the flow, i.e. it does not depend on the viscosity 
or the energy dissipation mechanism. The common entrainment assumption is that the 
entrainment velocity Ue is proportional to the typical velocity u inside the turbulent 
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zone (Morton et al. 1956; Turner 1986), usually the centreline velocity. The entrainment 
coefficient a — Ue/u is typically 0(0. 1), but the value is far from universal and depends 
on the choice of the typical length scale 5, the assumed shape of the velocity profile and 
can also depend on the initial conditions (e.g. Redford et al. 2012). 

Conversely, Corrsin emphasised a local perspective and suggested that the turbulence 
boundary is demarcated by a very thin viscosity-dominated laminar supcrlayer, whose 
local propagation velocity w„ towards the non-turbulent region is determined by two pa- 
rameters: the kinematic viscosity v and the rate of dissipation of kinetic energy e (Corrsin 
& Kistler 1955). Consequently, w„ oc u^ where u^ is the Kolmogorov velocity scale. The 
ratio of local entrainment velocity Vn to global (bulk) entrainment velocity Ug is given 

by 

!^oc^=Rc-l/^ (1.1) 

Me U 

and the Reynolds number dependence begs the question in which way the two views - local 
and global - are consistent. The dependence on Re seems to imply that both surface area 
and viscous diffusion adjust to the imposed global rate such that the small scale details of 
how the vorticity is transferred are somehow forgotten across interactions of eddies with 
a large hierarchy of sizes (Townsend 1976). By denoting the integral entrainment flux as 
Qp, the global perspective suggests that Qp — u^A, where A is the surface area based 
on the average distance of the turbulence interface to the core of the turbulent zone. 
From the local perspective Qe = VnS where S is the total surface area of the contorted 
turbulence boundary. Equating the two expressions for Qp results in Up/vn ~ S/A and 
therefore 

4 oc Re'/\ (1.2) 

A ' ' 

This means that S must be large to compensate a slow viscous transfer of vorticity and 
to cancel out the viscosity dependence (e.g. Tritton 1988; Sreenivasan et al. 1989). 

Probably the simplest setting where turbulence propagates into non-turbulent fluid is 
the case without any mean flow, which can be realized via oscillating grid experiments, 
e.g. Holzner et al. (2007, 2008); Holzner & Luethi (2011). The results obtained in such a 
flow showed evidence for the presence of a laminar superlayer at the boundary of turbulent 
flow regions. In particular, the analysis supported that S is indeed given by a strongly 
convoluted surface and accounts for a large entrainment flux with a small characteristic 
velocity comparable to the Kolmogorov velocity (Holzner & Luethi 2011). A similar 
picture, i.e. w„ ^ m,;, recently emerged from the experiments in a round jet of Wolf 
et al. (2012). In all the experiments and simulations, the probability density functions 
(PDFs) of the entrainment velocity indicated that there is a large distribution in the 
actual entrainment velocity. In this context the term laminar superlayer is unfortunate 
as it suggests that the flow is layered without notable fluctuations. Therefore, this layer 
will be termed the viscous superlayer (VSL) in the remainder of the paper. 

A somewhat different view emerges from direct numerical simulations of plane jets 
(Da Silva & Pereira 2008; Da Silva & Taveira 2010; Taveira & Da Silva 2013), a plane 
wake (Bisset et al. 2002) and experiments in a round jet (Westerweel et al. 2005, 2009), 
which focused on properties of the turbulent/nonturbulent interface (TNTI) . The TNTI 
seemed to be thicker than the VSL predicted by Corrsin; that is, the thickness of the 
TNTI was found comparable to the Taylor length scale A, rather than than the Kol- 
mogorov length scale 77, although it was conjectured that the VSL may be contained 
inside the TNTI region Bisset et al. (2002). Also, the physical picture promoted differs 
from the one proposed by Corrsin in that it suggests that turbulence propagates mostly 
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via transmission of turbulent (i.e. Reynolds type) shear stresses (Westerweel et al. 2005, 
2009), whereas it is the action of viscous shear forces in Corrsin's theory. 

One possible explanation for these differences in layer properties is the method by 
which the interface between the nonturbulent and turbulent fluid is identified. The inter- 
face is usually obtained by applying a threshold to a scalar field such as cnstrophy (Bisset 
et al. 2002; Holzner et al. 2007; Da Silva & Pereira 2008) or a high Schmidt number 
dye (Westerweel et al. 2005). By construction, this interface is artificial because the tran- 
sition between turbulent and nonturbulent fluid must occur smoothly over a flnite region. 
The assumption is that the results are insensitive to the precise value of the threshold, 
and papers tend to provide anecdotal evidence that this is the case. For example, the 
fact that the interface is sharp is commonly used as an argument for level independence 
of the qualitative findings, as long as the thresholds lies within an interval where the gra- 
dient is sharp e.g. (Bisset et al. 2002). In addition, in experiments and even in numerical 
data sets it is often difficult to vary thresholds over a span of several decades because of 
experimental (e.g. Westerweel et al. 2005; Holzner et al. 2007) or numerical (e.g. Bisset 
et al. 2002) noise. 

The aim of this paper is to perform a systematic study of the effect on the threshold 
value on the entrainment velocity and related statistics. Our analysis supports the exis- 
tence of a VSL over a large range of thresholds (^ 20 decades), with a smooth transition 
to the turbulent core region. The study provides insight into how the local and global 
turbulent entrainment are connected to each other and what it implies for the geometry 
of the bounding surface. We find that the VSL is a nearly one-dimensional layer with low 
surface curvature and both its area and viscous transport velocity adjust to the imposed 
rate so that the entrainment fiux is independent of threshold. We perform simulations 
of a temporal jet at Re = 5000 which has the advantage that it is well-documented 
(Da Silva & Pereira 2008), is a generic shear fiow and is attractive from the viewpoint 
that this fiow has two homogeneous directions. The present paper starts out with a the- 
oretical framework that describe the properties of the local entrainment velocity from 
a local and a integral perspective. Predictions are made based on a simple model and 
tested against the simulation data. Next, properties of the geometry and surface area 
of the VSL are presented. In the final section the results are brought in context with 
previous work and possible non-universality issues are discussed. In a companion paper 
the infiuence of mean shear and Reynolds number are analysed. 



2. Theory 

This section revisits the determination of the local entrainment velocity based on 
enstrophy budgets from a local perspective and complements it with an integral approach 
to the problem. Thereafter a simplified model is set out with predictions for the enstrophy 
transport across the VSL. 

2.1. Local entrainment velocity: classical approach 

We separate between turbulent and non-turbulent fiow regions by using a threshold on 
the enstrophy lo"^ = uJiUJi, where w^ is a component of vorticity. This defines a bounding 
surface separating the two regions. In a Lagrangian frame moving with the iso-enstrophy 
surface the isolevel will by definition remain constant and this property can be used to 
derive an expression for the propagation velocity (Holzner & Luethi 2011). We write the 
velocity of an isosurface element, tJ, as a sum of fluid velocity, it, and velocity of the area 
element relative to the fluid, V , that \s, v — u + V . The total change of w^ in a frame 
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of reference moving with an enstrophy isosurface element is then given by 

— . — + ..v.^^^ + y.v.^^o, (2.1) 

where the lower case A/ At is the total derivative following a surface element, and the upper 
case D/Dt is the material derivative which follows a fluid element. By defining a surface 
normal as n = Va;^/|Va;^| and the normal relative velocity component Vn = V ■ n, we 
obtain 

Substituting the enstrophy balance equation 

— f y j = i^V^ ( Y ) + '^'^J'^y ~ ^^^^ ■ ^^^ (2-3) 

into (2.2) and averaging over the isosurface {■) s, we obtain an expression for the average 
entrainment velocity w„: 



where 



t'i, (2.4) 



2ujiUjjSij\ Ti / i^ V cj^ \ r / 2iy'VoJi ■ 'Vuji 



Using the definition of n, the viscous term can be decomposed into a contribution due to 
curvature and normal transport through the following identity (Holzner & Luethi 2011): 

Vw2 = |Vu;2|V-n + n- V|Vu;2| (2.5) 

This identity will be used to quantify the role of curvature in §4.2. 

2.2. Local entrainment velocity: integral approach 

An alternative to the local approach described in § 2.1 is to integrate the enstrophy 
equation (2.3) over a time-dependent domain D{t) which has boundary velocity v and 
use the Reynolds transport theorem, resulting in 

D Jd 

As the surface normal h points into the turbulent region, the appropriate volume under 
consideration comprises the irrotational region. We can formalise this by defining a con- 
trol volume D = H{1 — cj^/cjq) where H is the Heaviside function and ujq is a enstrophy 
threshold. As w^ is then by definition constant on the surface dD, the equation above 
simplifies to 

V-ndS^^(^[ ^dV-4 v(^]-ndS 
9o oji \dtJij 2 Jan \2 J ^2 6) 

— / ujiUjjSijdV + V I Vwi • VcjidV" 
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Definingf 

dV 



V • ndS = VnS, S= (f> d5, V = dV, (2.7) 

dc JdD JdD Jd 

and introducing an average over the volume (•)y, (2.6) can be rewritten as 



{^^)v\ _ 2i./dc.2\ 2V ( d {io^)y 



4\dn / g Suj^ 



Now, because D spans the entire nonturbulcnt region, it is expected that (lj'^') /cjq <C 1, 
and therefore 

Vn'^^vf +vj + vf +v^, (2.8) 

where 

r_ 2V d {iV^)y 

^' Suj^dt 2 ' 
^ 2V 



v? = 


2u /duj^\ 
^0 \ dn Z^' 


<- 





^c.2 



2.3. ^ model for enstrophy transport in the viscous superlayer 

In the VSL, the evolution of enstrophy is governed by molecular processes (Corrsin & 
Kistler 1955; Holzner & Luethi 2011), i.e. \v'^ /vn\ ^ 1- Assuming that the local curvature 
is small on average and multiplying by dcj^/dn, Eq. (2.4) then becomes 

da;2 d^w d / du;\ , , 

-— w„ + 2vu:-—r = 2uj— VnijJ + v-— = 0. (2.9) 

dn dn^ dn \ dn J 

Integrating this expression and using that at n = — oo, both w = and duj/dn = 0, we 
obtain 

w„w + t/^=0. (2.10) 

dn 

Assuming that w„ is constant in the VSL, the square of the solution to (2.10) is 

— = exp , 2.11) 

ujI \ V J 

where w^ = u]{n.r) is a reference value for enstrophy. Hence, the enstrophy is expected to 
drop off exponentially in the VSL, provided that w„ is constant (see §4.2). 

The model solution (2.11) can be used to predict the magnitude of the enstrophy 
transport terms. For the local approach we expect 

dn^ \ dn ) 

f / du}\ / diu^ ^ 

vi^2vi --] -— 
\dn J \ dn 

f The first equation holds because £ u ■ ndS = for an incompressible fluid. 
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Entirely consistently, we expect for the integral approach that 



.."D 



= 2w„, 



uj^ dn 



3. Simulations 

The start situation for a temporal plane jet is a fluid layer that is quiescent except 
for a thin region — feo < y < bo 'm which the streamwise velocity u is nonzero, and is 
homogeneous in the two other directions x and z. Here, bo is the initial jet width. It 
follows from continuity that the volume flux q = J udy remains constant for this flow 
throughout the jet's transition to turbulence and subsequent growth due to turbulent 
entrainmcnt. Assuming that the Reynolds number Re 3> 1, the only relevant parameters 
are the initial volume flux qo and time t which suggests self-similar behaviour, with the 
jet halfwidth b and centreline velocity ii scaling as 6 oc ^/qot and ii oc \/qo/t, respectively. 

The simulation domain is a cuboid of size 245o x 366o x 246o, which is three times larger 
in all directions than the domain used in Da Silva & Pereira (2008). The larger domain 
facilitates much longer simulations, thereby allowing not only the first moments but 
also the turbulence to reach an equilibrium and has the added advantage of improved 
statistics because of the larger area spanned by the two homogeneous directions. The 
simulation considered here is for Re = 2qo/v = 5000. The resolution of the simulation 
is 1024 X 1536 x 1024 which is sufficient to ensure that all active length scales in the 
turbulence are fully resolved. We define a reference time-scale t* = b^/qo and simulate 
for soot*. All statistics before t/t* ==150 are discarded to ensure that the turbulence has 
time to reach a dynamic equilibrium. 

Following Da Silva & Metais (2002); Da Silva & Pereira (2008), we use the initial 
condition 



"(y,0) = ^ 



l + tanh'^°-"^" 



(3.1) 



200 

where C/g is chosen such that J udy = qo- The initial condition is then seeded with uni- 
form random noise of such amplitude that the enstrophy level due to noise is 8 percent 
of the maximum average enstrophy level [which is (J7o/(46'o))^]. The code used for di- 
rect numerical simulation is based on fully conservative second order finite difference 
operators in space (Verstappen & Veldman 2003) and uses an adaptive second-order 
Adams-Bashforth time integration scheme. More details can be found in van Reeuwijk 
et al. (2008). 

As the statistics shown in the next section depend heavily on budgets for the enstrophy, 
special care is taken to ensure that the budgets are calculated consistently with the 
numerics used. To achieve this, a mimetic (Hyman & Shashkov 1997; van Reeuwijk 2011) 
curl operator is defined such that it satisfies the identity V x Vp = up to machine 
precision, where p denotes pressure. In order to ensure calculation fully compatible with 
the numerical method, we do not manually discretise (2.3), but instead make use of the 
following identities 

uJiUjjS^j = (m • V)— a; • (V X {u ■ V)m), (3.2) 

i^ViJi ■ Vlj, = i/V^ ( ^ j - i/w • (V X V^tt), (3.3) 



Jpf pnf.rnin.m.pni- from th.p -ii-isipmi.^ fiiinprlaiipr tn thp tiirhiilpmt pnrp rptrinr 






(a) 



■y = 



Cn; 




Figure 1. (a) Aa; normalised by the Kolmogorov scale r\\ (b) Re\ against time. 

which are then also enforced on the discrete leveL Taking (3.3) as an example, one can 
calculate the first term on the right hand side directly by using the routine for scalar 
diffusion on a;^/2; the second term can be calculated by taking the discrete curl of the 
viscous term in the momentum equation, and then taking the scalar product of the result 
with the vorticity components. 

As the temporal jet has a nonzero mean velocity in the cc-direction, it is important to 
ensure that the identity ^ u ■ fidS = is also satisfied on the discrete level. Indeed, this 
identity was used explicitly to derive (2.7). This can only be achieved if the thresholding 
algorithm identifies entire cells to be either inside or outside the turbulent region. Indeed, 
we found that if we used bilinear interpolation to construct an isosurface - which is in 
principle a better representation - the various interpolations required could lead to very 
significant deviations in the calculated entrainment velocity. 



4. Results 

4.1. Bulk flow properties 

Fig. 1(a) demonstrates that the grid resolution is appropriate for the problem under 
consideration. Shown is the grid spacing normalised by the Kolmogorov length scale 77 = 
(i^^/e)^/"' based on the centreline dissipation rate e{t) = e{y = 0, t). The overbar denotes 
averaging over the two homogeneous directions and over 10 t*. The dissipation rate has 
its maximum at y = and dissipation rates will be much lower at the jet boundary, which 
implies that the simulation is even better resolved there (dashed and dash-dotted lines). 
As can be seen, the simulation becomes better resolved in time because 77 oc yt, which 
can be inferred by using that e ex u'^/b, as is confirmed in Fig. 3(d). Fig. 1(b) shows the 
evolution of the Taylor Reynolds number, defined as Re\ = (20/3)e^/v^ (e.g. Tennekes 
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Figure 2. Energy density (left) and dissipation spectra (right) at y = 0. 



& Lumley 1972), where e{t) = e{y — 0, t) is the turbulent kinetic energy. As judged from 
Re\, the turbulence reaches an equilibrium value of about Re\ w 100 after t/t* — 50. 

Energy density spectra for the plane y = 0, averaged over shells of wave number 
a//c^ + k1 and a time interval of lOi* are shown in Figure 2. The spectra indicate that 
there is a range of active scales and that there is a separation of scales as is evident 
from the formation of a k~^'^ spectrum (left) and the peak in the dissipation spectrum 
(right). Sixteen spectra are shown for t/t* > 150 and the collapse clearly demonstrates 
the self-similarity of the flow under consideration; even though the spectra change in 
time, the normalisation with 77 and e cause a full collapse of the data. 

As the velocity profile is symmetric around y ~ 0, the jet half-width b was inferred 
from the average of the values of b for which u{b, t) = u{t)/2 and u{—b, t) = u{t)/2. For 
all profiles shown, use has been made of the symmetry (or anti-symmetry) in the profile 
to further improve the statistical accuracy. Shown in Fig. 3(a) is the scaling of &^ with 
time, which as expected from dimensional arguments is linear; the red dashed line is a 
linear fit. The normalised mean velocity m and momentum flux v'u' are shown in Fig 
3(b,c). These profiles were scaled and then further averaged over four contiguous time- 
intervals (an effective average over 40 t*). The profiles are convincingly self-similar. Self- 
similarity of turbulent quantities is demonstrated in Figs 3(d-f). As discussed earlier, the 
balance between turbulence production and dissipation suggests that e oc u^ /b (Fig. 3(f)), 
which in turn implies that rj cc b (Fig. 3(d)). The profiles for the turbulence kinetic 
energy (TKE) e and the dissipation rate e (Figs 3(e, f)) show once more a reasonably 
good collapse, with e peaking at y/b ~ 1 (where shear-production is maximal) and the 
dissipation rate e peaking at the centreline. Consistent with recent work by Redford et al. 
(2012), we observed quite substantial variations in the growth rate of b between different 
simulations despite a convincing self-similar behaviour in all of them. This points to 
non-universal self-similar behaviour (Redford et al. 2012). 

4.2. Entrainment velocity and budgets 

A cross-section of the enstrophy field at constant z is shown in Fig. 4, together with a 
number of enstrophy isolines. Bisset et al. (2002); Da Silva & Pereira (2008) use a thresh- 
old value Wo = 0.7L/o/(2foo) which corresponds with loq/uj^ « 0.1, where we have defined a 
reference enstrophy level w^ = U^/bQ. At this enstrophy level, the turbulent /nonturbulent 
interface is highly contorted and has 'holes'. Also shown are two enstrophy thresholds 
that are much lower, namely Wg/w^ = 10""'^^ and 10~^. What is striking is that these 
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Figure 3. Self-similarity of the temporal plane jet. (a) dependence of h^ on time; (b) average 
velocity u\ (c) turbulent momentum flux v'u'; (d) dependence of h/rj on time, (e) kinetic energy 
e; and (f ) the dissipation rate e. 
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Figure 4. A field in cross-section showing ^" \og{u)'^ / u'r) at t/t* = 150. Also shown are various 

isocontours of enstrophy. 



low enstrophy levels are in close proximity to the high enstrophy regions, indicating a 
very quick drop-off of enstrophy levels near the jet edge. At these low threshold values, 
there are no more holes visiblef in the figure, and although the surface remains contorted 
because of the large scale vortices distorting the flow, the enstrophy isosurfaces appear 
to form nearly one-dimensional layers with relatively small curvature. 

The instantaneous budgets of enstrophy were calculated for the entire 3-D field every 
5 t* and then used to calculate the terms in (2.4), (2.8) for 37 thresholds in the range 

f Note that what may appear as holes on the Fig. (i.e. in 2D) are fluid portions that are 
connected in 3D to the outer irrotational region. 
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Wq/w^ € [lO"^**, 10"]. Simultaneously, the volume V was recorded for each of the threshold 
values, which enabled an independent calculation of the propagation velocity w„ using 
(2.7). Shown is the data for t/t* > 150. 

The terms comprising the local entrainment velocity v^ (2.4) normalised by the directly 
measured volume-based entrainment velocity w„ are presented in Fig. 5(a). For loq/uj'^ < 
10^"'', the predicted propagation velocity u^ (squares) excellently matches the actual 
propagation velocity. The poor predictions for cuq/lu'^ > 10~^ are related to an insufficient 
temporal sampling frequency creating large errors in the calculation of Vn ~ 5V/St (Eq. 
2.7) at high enstrophy levels. 

Consistent with earlier findings, the inviscid contribution v^ (circles) does not play a 
role in the VSL (Holzner & Luethi 2011). For oJq/lu'^ > 10^'"', the inviscid terms become 
dominant very rapidly, and note that the enstrophy level employed by Da Silva & Pereira 
(2008) and others is cjq/cj^ — 10^^. Also shown in Fig. 5(a) are the theoretical predic- 
tions from §2.3, namely v^ — 2w„ and w£ = — w„ (both displayed by dotted lines). The 
predictions are in good agreement, although the observed magnitudes of v^ (downward 
triangles) and w„ (upward triangles) are a bit larger. The plots of vl reveal a slight 
decrease at lower threshold levels, indicating that there is a small systematic error in the 
calculation of v^. 

Fig. 5(b) shows an identical calculation but using the integral approach outlined in 
§2.2. As for Vl, the propagation velocity vj is in excellent agreement with w„ for luq < lo^. 
The temporal and inviscid contributions, vj (diamonds) and uj' (circles) respectively, 
are negligible and the propagation of the enstrophy isosurface in the VSL is caused by 
viscous effects. It is noted that for this indicator, the contributions vY and wf seem to 
become independent of the threshold level below clIq/cj^ < 10^'"'. Furthermore, vi/vn 
becomes constant for small threshold values, although the ratio is slightly smaller than 
1. 

One of the main assumptions made in the derivation of the theoretical model was 
that the curvature of the isosurface is small. The reasonably good agreement with the 
theoretical model in Fig. 4 supports this assumption but using Eq. 2.5 it can be validated 
explicitly. Figs 5(c,d) show the contribution to v'^ by curvature (red crosses) and diffusive 
transport in the direction of the surface normal (blue triangles). For both vl as u/, the 
curvature term becomes negligible for uj'^lujf, < 10^^. This provides confirmation that 
the theoretical model in §2.3 is a reasonable description of the processes governing the 
VSL. 

Having established that both estimates of the local entrainment velocity vl and vj are 
in good agreement with the actual entrainment velocity u„, we study the dependence 
of the entrainment velocity on the threshold value Wq. The bulk entrainment velocity is 
defined as Ue = L~^dV/dt where L — 24feo, and is therefore related to «„ (Eq. 2.7) by 
Ue = —{S/L'^)vn- In Fig. 6 we show the dependence of «„, vl and vj on ujq. 

There is a clear dependence of Vn on the enstrophy threshold when normalising by 
Ujj [Fig. 6(a)]; isosurfaces for the very low enstrophy thresholds propagate faster than 
those with higher thresholds. Indeed, vl is nearly twice as high at Wq/o;^ = 10^^** than 
at ujq/uj'^ = 10^^. Hence, although w„ is of the same order of magnitude as Ujj, the de- 
pendence on cjq suggests that it is not merely the viscosity v and dissipation rate e that 
determine the propagation velocity in the VSL. Another striking feature is that w„ be- 
comes zero around Wq/w^ « 10~^ for the timespan under consideration and negative for 
uig/uir > 10~^. Hence, high enstrophy regions move towards the jet centre and the low 
enstrophy regions move outwards. This can be explained by using the relation between 
enstrophy and the dissipation rate which is valid for isotropic and homogeneous turbu- 
lence e — vuj'uj'. Using the self-similarity of e it follows that w'w' « u^ /{vb) ex t'"^ . Hence, 
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Figure 5. the local entrainment velocity decomposed into various contributions and normalised 
by the actual propagation velocity «„. v^ [A], v^ [o], v^ [V], u^ [0], total [D]. (a,c) vl\ (b,d) 
vi. Figs (c,d) show the decomposition of v^ in a contribution due to curvature (red crosses) and 
normal transport (blue triangles). 



if one would assume that ijj'uj' has a self-similar profile and pick a reference threshold 
Wg that is followed in time, it would be seen to move towards the jet core. This applies 
to enstrophy levels where turbulence is developed and u„ is positive. Towards the VSL, 
at low aj2 levels, viscous transport is diffusing enstrophy outwards and u„ is negative. 
Note that the gradient of w„ over uP in Fig. 6(a) is always negative meaning that the 
enstrophy profile is flattening over time. 

Fig. 6(b) shows the bulk entrainment velocity Ue normalised by the centreline velocity 
amplitude u. The standard entrainment assumption (e.g. Turner 1986) is Ue — ait and the 
ratio plotted is the entrainment coefficient a. As can be seen, a is much less dependent 
on Wg < 1^2 than Vn/uri- If a were entirely independent of ojq, it would imply that the 
entrainment flux Qe is constant; indeed, by using that Ue = v„S/L'^ = Qe/L^ it would 
immediately follow that Qe is constant if Ue is independent of ujq. Even though a small 
dependence on a is discernible, this result does indicate that to first order S (x v~^, and 
only if Vn oc m^ is independent of Wg do we expect that S oc Re ' , independently of Wg, 
i.e. the classical view advocated by Corrsin. 

Full independence of the entrainment velocity Ue of the enstrophy threshold ujq can 
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Figure 6. Entrainment velocities as a function of threshold, (a) u„ normalised by it,,, (b) u^ 
normalised by u. (c) normalised entrainment coefficient qq- 



be achieved by taking into account the position of the interface. Selfsimilarity implies 
that u — uf{£,), where / is the universal velocity profile and ^ — y/b the selfsimilarity 
variable. The entrainment assumption is Ueo = ctQU, where Weo is defined as Ueo = d5/dt. 
As mentioned in §1, other definitions of b will result in different values for a. This is 
straightforward to see by fixing ^ at a value differing from unity. Indeed, for an alternative 
width y = ^b, we obtain 



ao£,u, 



(4.1) 



dy .d6 
"^=dt=^d^ 

indicating that the effective entrainment rate is u^/u = fap. Hence, if there is a depen- 
dence of the average interface position {y)g on the enstrophy isosurface this will create 
a trend in Ue- In Fig. 6(c), we have plotted ao = Ue/{£^u), where ^ = {y)g/b. As can 
be seen, both u„ as uj are constant as a function of the threshold value. There is a 
downward trend for ul which is consistent with Figs 5(b) and 6(b), but as u„ is directly 
calculated from V , we consider this data to be more trustworthy. 

There are a two distinguishing enstrophy threshold values to consider at any time, 
namely (i) the enstrophy level at which w„ = and (ii) the enstrophy level for which en- 
strophy production becomes negligible, diagnosed by the operational criterion \v'^ /vn\ < 
0.05. The former demarcates a threshold which separates expanding areas (w„ < 0) from 
shrinking areas (w„ > 0) and the latter demarcates the start of the VSL. These threshold 
levels are plotted as a function of time in Fig. 7 for vl (sohd lines) and vj (dashed lines). 
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Figure 7. A map of the different regions as identified by the thresholds on enstrophy; vl 

(solid lines) and vi (dashed lines). 



Although the exact values obtained by the two estimates differ, the trends are very con- 
sistent. As time progresses, the threshold level for which v„ = moves to lower and lower 
levels, in accordance with the decay expected from self-similarity The threshold level de- 
marcating the beginning of the VSL remains approximately constant at ujq/uj'^ « 10~^ 
for t/t* > 150. 

Using these two threshold levels we can distinguish three regions: (i) the core region, 
where inviscid production, viscous transport and dissipation are all important; (ii) the 
intermediate region, which can be compared to the buffer layer of a wall-bounded turbu- 
lent flow in the sense that both viscous and inertial effects are important; and (iii) the 
VSL, which is dominated by viscous transport and dissipation. The results presented in 
Da Silva & Pereira (2008); Da Silva & Taveira (2010); Taveira & Da Silva (2013) are 
all obtained with thresholds values corresponding to the core region. Note that, similar 
to the viscous sublayer in a wall-bounded flow, the VSL can neither be classified as tur- 
bulent because viscous effects are dominant nor as laminar because there are significant 
fluctuations in the layer due to external influences. Fig. 8 demonstrates how close the 
turbulent core region is to the VSL. The core region occupies most of the turbulent fluid 
and contains many "holes", whereas the VSL forms a thin strip in the white region im- 
mediately adjacent to the isoline marked as |w^/wn| < 0.05. The intermediate region is 
very thin, and seems to be thinner ahead of the turbulent events travelling from left to 
right on the figure, and thicker in their wakes. 

4.3. Geometry 

The dependence of volume V and surface area S on the threshold value ujq can be used 
to obtain information about the average distance from one enstrophy isosurface to the 
next, thereby getting an impression of the distance between different regions of the flow. 
The average distance n, the volume Vr for which w^ > Wgt and the surface area S are 



t The sum of V and Vr is exactly half of the simulation domain volume. 
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Figure 8. A field in cross-section showing ^ log(i:i; /cj^) at t/t* — 150. Clearly visible is the 
close proximity of the core region to the VSL. 



related by 



Introducing AVr — Vr-i+i — Vt-, 



dVr 
dn 
and An = 



= S. (4.2) 

rii+i — Hi as the difference in volume and 



distance respectively between two adjacent enstrophy thresholds Wq.j and (jJQ.i+i, the 
average distance can be approximated by An sa 2AVT/{Si + Si+i). We set n = at the 
start of the VSL at wg/w^ = 10~5. 

The dependence of n on ljq is shown for t/t* > 150 in Fig. 9(a). The collapse of the 
profiles for different times shows that the result is quite robust. It is striking how close 
essentially irrotational regions are located to regions with very high enstrophy levels on 
average. Indeed, the core region is on average about IO77 away from the VSL. Moreover, 
the full 24 decades in enstrophy levels are separated by 25?]. Also shown is the analytical 
solution (2.11) (dotted line), expressed as Wq/o;^ = exp(Q;^(n — nr)/r]) with a,, = 3.7. It 
is evident that (2.11) is a reasonable approximation for coq/oj^. <C 1, which explains the 
agreement of the simulation data with the predictions for u^ and v^ . However, as can be 
seen, the shape of the n profile deviates from a straight line in a semilog plot because of 
the influence of surface area S* on w„ . Indeed, the model assumes that w„ is constant and 
this has clearly shown to be not the case for this flow. Fig. 9(b) shows the normalized 
distance for a number of thresholds as a function time. After t/t* = 150, the distances 
become nearly constant, indicating that the jet edge has reached a dynamic equilibrium 
and the distances between isosurfaces scale with 77. 



5. Concluding remarks 

The results presented in this paper show convincing evidence for the existence of the 
VSL. Consistent with earlier work, we find that inertial processes are negligible in the 
VSL and that this layer is discernible for nearly twenty orders of magnitude in enstrophy 
threshold, a range corresponding to an average distance of 15 77. The VSL was present for 
oJq/uj^ < 10^^ and was practically independent of time for the duration of the simulation. 
The simple theoretical model derived in §2.3 is in good agreement with the data and 
shows that the contribution of the viscous transport term amount to 2w„ and the viscous 
destruction term to — u„. 

Whilst Vn was of the same order of magnitude as u,,, there was a clear dependence 
of the value of Vn on the enstrophy threshold. This suggests that Corrsin's dimensional 



Jet entrainment: from the viscous superlaver to the turbulent core region 



15 



e 




e 



10-24 10-18 10-12 10-6 10° 

,2 /, ,2 



^oK 





ullu^l = lQ-]_ 



ullul = 10- 



100 200 

t/t* 



300 



Figure 9. (a) thickness as a function of threshold value; b) thickness as a function of time. 



arguments may not convey the entire story of the outer fringes of a turbulent zone. Indeed, 
we found that the classical inviscid entrainment assumption Ue — aii best described the 
entrainment in the VSL as Qe — v^S is practically independent of the threshold value. 
This implies that both viscous transport velocity and surface area adjust to the imposed 
global rate so that the outcome is independent of threshold. As the VSL becomes less 
contorted when moving out of the turbulent region w„ needs to become larger to maintain 
a constant entrainment flux. Note that the inviscid entrainment assumption works for 
the global behaviour, but the local physics of the VSL are all but inviscid and governed 
by viscous diffusion of enstrophy. 

The small dependence on Wg that was observed when plotting Ue/u could be explained 
by taking into account the average position {y)g /b of the interface. Using that the VSL 
spans ~ 15r7 and b ~ 150r/ the expected difference is about 10% for the case considered 
here, which is consistent with Fig. 6(b). Employing Corrsin's assumption that the thick- 
ness of the VSL is ~ ?7 and noting that b ~ r]Re ' , it follows that as Re tends to infinity, 
the ratio {y)g /b — > will become constant throughout the VSL. Hence, this appears to be 
a low-Reynolds number correction. 

Three regions were observed for the flow under consideration. For roughly Wq/w^ < 
10~^, we observed the VSL which is characterised by a propagation velocity only de- 
pending on viscous processes. The core region was categorised as that region of the flow 
for which u„ > 0, i.e. regions that become smaller in time and that are responsible for 
the uj'uj' (X i-^ behaviour. The enstrophy threshold at which this region was observed 
was a decreasing function of time. In between the VSL and the core region there was an 
intermediate region providing a smooth transition from one region to the other. 

This study suggests that the reported differences in the thickness of the TNTI and 
the VSL are due to different choices of the threshold. Indeed, Da Silva & Pereira (2008); 
Da Silva & Taveira (2010); Taveira & Da Silva (2013) use a threshold value which is 
associated with the core region. The distance from the core region to the VSL is about 
10?7 for the Re under consideration, which is a value not atypical to that mentioned in 
Da Silva & Pereira (2008). The fact that the core region is fully turbulent may provide 
an explanation for why the layer thickness is found to scale with the Taylor length 
scale A. We should point out that another possible explanation for the difference might 
be that Da Silva & Pereira (2008) use conditional statistics that are dependent on the 
inhomogeneous direction only - and not along isosurfaces as was done in this paper. 
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Further work, particularly a study of the Re dependence is necessary to be able to make 
a firmer statement about the difference in boundary layer scaling. 

Finally, it may be argued that the term turbulent-nonturbulent "interface" is slightly 
misleading, as it suggests (i) a surface separating two regions; and (ii) a sharp distinction 
between the character of the flow in these two regions. Consistent with earlier work, we 
find that the transition from turbulent to non-turbulent occurs over a thin but finite 
region, in which the flow changes smoothly from one region to another. Any particular 
threshold value to demarcate the turbulent/nonturbulent interface is then by definition 
construed. Furthermore, the current work suggests that the TNTI comprises both the 
intermediate region and the VSL, confirming the conjecture that the VSL is contained 
inside the TNTI region made by Bisset et al. (2002). 
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